A spatially explicit model for tropical tree diversity patterns 
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I. INTRODUCTION 

On long time scales — evolutionary time — the particu- 
lar characteristics of species evolve widely from birth to 
extinction, while throughout ecological time most species 
persist sometimes for decades or even centuries [l], 
Their organization in complex networks of local com- 
munities and extended metacommunities is subject to 
abundant variations in terms of relative distributions. 
Usually, species maintain themselves with constant abun- 
dance but, from time-to-time, processes such as invasion 
or succession produce fast regrouping throughout a par- 
ticular site or region 0|. 

Macroecology 0, [3, IH studies the relationships be- 
tween organisms and their environment at large spatial 
scales in order to characterize and explain universal sta- 
tisticalpatterns of abundance, distribution and diversity 
@, 0, Q- A variety of models and methods from sta- 
tistical physics are appropriate for the study of ecosys- 
tem dynamics, in which the main "entities" are either 
many individual organisms within populations, or many 
species within local, regional or continental communities 

BE, El, EEl. 

The remarkable regularities in patterns of data on how 



species originates, persist, assemble in groups, and even- 
tually go extinct, suggest the existence of general mech- 
anisms from which the biodiversity and the structure of 
ecological communities originate 0, 0, 0] ■ The difficul- 
ties in explaining the biological diversity of such systems 
originate from the very different spatial and temporal 
scales, starting with the evolution and biogeographic dis- 
tribution of species and ending with individual births and 
deaths in local communities. Present approaches for de- 
scribing the dynamics of population genetics and ecology 
consider either the importance of the genetic fitness or 
the influence of random events governing birth, death and 
migration phenomena p[. There are thus two main, but 
conflicting perspectives on the nature of ecological com- 
munities: the niche-assembly and the dispersal-assembly 
perspective [3, 0| ■ In the present work we intend to use 
this second approach, constructing a simple spatially ex- 
plicit stochastic model for describing the relative species 
abundance and the spatial diversity patterns of trees in 
a tropical forest [13, EE EH]- 
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II. NEUTRAL MODELS 

Nowadays scientist argue that ecological systems of 
similar species that compete with each-other only for the 
limited amount of natural resources can be successfully 
approached by means of neutral models. One of the first 
and most widely used neutral models is that of Stephen 
Hubbell [3, [T3|, who assumed in his theory that within 
groups of ecologically similar species, individuals fill the 
landscapes up to a point of saturation and the dynamics 
of all these species is governed by the same birth-, death-, 
migration- and mutation-rates. This simple neutrality 
principle allows us to study a wide variety of systems by 
means of stochastic computer simulations and to use the 
methods of statistical mechanics to get analytical results 

[H, eh mi, m. 

Tropical forests [3, [H, M, HI contain a wider variety 
of tree species compared to northern coniferous forests (2 
orders of magnitude larger) , and it is believed that such 
systems are much closer to the neutrality assumption. 
Since the number of mature individual trees in equal- 
area samples of each forest type is almost the same, the 
increase in the abundance of a particular species needs to 
be compensated by a decrease in another species's pop- 
ulation. 

In order to develop an appropriate model, one needs 
to take into account the whole set of factors influenc- 
ing community composition: births, deaths, immigration 
and (on longer time scales) speciation. These processes 
interact in a complex manner and produce the commonly 
observed empirical patterns of diversity and distribution 
of abundances. The model should prescribe, for each gen- 
eration, a method for choosing the species identity of the 
individuals that die and those that will occupy (by birth 
or immigration) the vacancies created. A simple method 
which describes fairly well the above mentioned processes 
is the one developed by Hubbell and it is called "zero-sum 
ecological drift" P,LL3]- This method implies that the re- 
placing species are drawn at random from the existing 
community of species. The "ecological drift" is connected 
to the random replacement process and assumes that all 
individuals, regardless of species, have equal probabili- 
ties of giving birth, dying, immigrating or acquiring a 
mutation to generate a speciation event. The ecological 
drift does not imply an equal chance for each species to 
fill a given vacancy. Obviously, the greater probability 
of being drawn into a vacancy belongs to more abun- 
dant species. A remarkable observation is that individu- 
als are equal, but species, as collective entities, are not. 
In Hubbell's model the ecological drift, with no ad- 

ditional mechanism is enough to produce the patterns of 
species abundance and diversity observed in nature. In 
particular, the community of tropical trees is very well 
described by these neutral models and neutrality will be 
also the basis of the approach considered in this work. 
Contrary to earlier modeling efforts [3, here we 

will continue the novel idea introduced by Zillio et al. 
[25[, and a spatially explicit model resembling the clas- 



sical voter model will be studied. In our approach we 
will consider, however, a more complex, two-parameter 
version, which is not suitable for an analytical study. We 
will argue why it is important to make the model con- 
sidered in [25] more complex, and we will determine the 
best parameters for optimally reproducing the measured 
macroecological patterns. 



III. RELEVANT MACROECOLOGICAL 
MEASURES 

In order to describe the statistics of species sizes in 
meta- or local ecological communities, the main mea- 
sure which is usually considered is the Relative Species 
Abundances (RSA) and the Rank- Abundance distribu- 
tion. Provided our data is spatially accurate, and we 
know the location of each individual, we can also inves- 
tigate the species-area scaling and the auto-correlation 
function for the individuals of a given species. We will 
discuss briefly these measures here. 



A. Relative species abundance distribution (RSA) 

The relative species abundance distribution is intro- 
duced for characterizing the frequency of species with 
a given abundance [H, Hf|. Three different types of 
plots are generally used (Fig. [1]) for representing species- 
abundance distribution. Historically the first, and thus 
also the most widespread, representation of species abun- 
dances is due to Preston [2fl|, who considered exponen- 
tially increasing abundance intervals (i.e. a 4 , i — 1, 2, . . . 
with fixed a), and plotted the number of species found 
within these intervals as a function of i. Preston's plot is 
motivated by the fact that abundances can vary in wide 
limits and there are relatively small numbers of abundant 
species. By considering fixed length abundance intervals 
one would get large statistical fluctuations on the tail of 
the curve. For widely different communities, the Preston 
plot has a gaussian-like shape. As an example for such a 
representation, on Fig. UK we illustrate the shape of RSA 
for the Barro Colorado tropical forest tree census. In this 
case however, the curve does not show a clear gaussian 
shape. 

A second way of representing the species abundance 
distribution is arranging the species in decreasing or- 
der by their abundances and plotting the species's rank 
versus its abundance (rank-abundance plot) on a log- 
log scale [Tol . [l2l [T3 |. This type of representation is 
inspired by several abundance studies in sociology and 
economics, leading to the very general Pareto-Zipf distri- 
bution. For the Barro Colorado tropical tree census the 
rank-abundance plot is illustrated on Fig.QJ). 

The third and mathematically most rigorous way of 
representing the species-abundance distribution is plot- 
ting the mathematical distribution function p(s), i.e. the 
probability density for finding a species with a given s 
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abundance. For most of the neutral-like communities 
this distribution function has a tilted J shape on a log-log 
scale. As an example, on Fig. [Tfc we plot this probability 
density for the Barro Colorado tropical tree census. It is 
worth mentioning that p(s) can be derived from the Pre- 
ston plot by dividing the number of individuals in each 
interval with the length of the interval, and plotting on 
log-log scale this quantity versus the mean abundance in 
the given interval. It is also easy to realize that the rank- 
abundance plot is related to the cumulative distribution 
function corresponding to the p(s). 
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FIG. 1: RSA for the 50 ha Barro Colorado Island tropical 
tree census. Three different ways of presenting RSA: (a) the 
Preston plot, (b) the rank-abundance plot, and (c) the proba- 
bility density function (the dashed line indicates a power-law 
with exponent —1). 



B. Species-area relationship (SAR) 

Generally, the number of detected species does not 
scale linearly by increasing the size of the sampled ter- 
ritory. In order to characterize this dependence, the 



species-area relationship is studied. This is done by con- 
sidering larger and larger territories, and counting the 
number of species present within these areas 0, UHH]. 
For a better statistic, sometimes an average species num- 
ber is calculated on several territories with similar areas. 
The dependence of the average species number on the 
size of the sampling area is usually a power-law, with ex- 
ponent in the range 0.15-0.5. As usual, the presence of 
the scaling is visually illustrated by using a log-log scale 
like in (Fig.EJ). 
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FIG. 2: Species-area scaling for different BCI census data. 
The continuous line suggests the power-law fit in the limit of 
large areas. Scaling exponents in this limit are in the range 
0.15-0.16 for the different census years. Both scales are loga- 
rithmic. 



C. The auto-correlation function 

In case the individuals are restricted to a given spa- 
tial position, like trees for example, one can characterize 
the spatial distribution of the individuals from a given 
species using an auto-correlation function. This function 
is defined as the usual auto-correlation function in sta- 
tistical physics, and characterizes the probability density 
that for a given individual an individual from the same 
species exists at a distance r. The way we construct the 
C(r) auto-correlation function is the following. First a 
check-board type uniform mesh is considered, and the 
territory is divided into small, square-like domains (la- 
beled by coordinates i, j) and the number of individuals 
from the considered species is determined in each domain 
(Nij). The auto-correlation function for a relative coor- 
dinate p, q (Cp^ q ) is calculated as 



((N itj - (JV))(2V, 



i+p,j+q 



<*))>, 



(1) 



Here (N) denotes the average number of individuals from 
the considered species in the constructed domains: (N) — 



Since there is no reason to assume that the 



distribution of the individuals is non-isotropic, we can 



calculate the average of the C Pj<2 values for all p, q values 
that are inside a ring with radius r and width Ar (r < 
\/p 2 + q 2 < r + Ar) , considering a reasonable small Ar 
value. 

C(r) = (C p , q {r < < r + Ar}) pq (2) 
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FIG. 4: Density of trees in the 1990 BCI census. The 
grayscale code indicates the density, darker squares indicate 
higher tree densities. 



FIG. 3: The C(r) spatial correlation function for the species 
that is most abundant among trees with a diameter larger 
than 10 cm (Trichilia tuberculata) . Both scales are logarith- 
mic. Figure a. shows the results considering all individuals 
and Figure b. is obtained for trees with diameter larger than 
10 cm. 



IV. EXPERIMENTAL DATA: THE BARRO 
COLORADO ISLAND DATASET 

The experimental results used for comparison are from 
a detailed 50 hectares tropical forest tree census in Barro 
Colorado Island (BCI), realized by the Smithsonian Trop- 
ical Research Institute, Center for Tropical Forest Science 
(CTFS) [li]. BCI is located in the Atlantic watershed of 
the Gatun Lake (Panama) and was declared a biological 
reserve in 1923. It has been administrated by the Smith- 
sonian Tropical Research Institute since 1946. From the 
viewpoint of ecological studies, this island is ideal, be- 
cause it is covered with a rain forest that is still unper- 
turbed by human beings. The flora and fauna of Barro 
Colorado Island have been studied extensively and in- 
ventories have reported 1369 plant species, 93 mammal 
species (including bats), 366 avian species (including mi- 
gratory), and 90 species of amphibians and reptiles. The 
tropical tree census was performed only on a small part 
of this huge island, precisely on a 1000 x 500 meters (50 
hectares) area. The first census was completed in 1982, 
revealing a total of approximately 240,000 stems of 303 
species of trees and shrubs more than 1 cm in diameter at 
breast height. This CTFS program has an important ad- 
vantage, because in each census, all free-standing woody 
stems at least 10 mm diameter at breast height are iden- 
tified, tagged, and mapped, and hence accurate statistics 
can be made. Data is publicly available for several years: 
1982, 1985, 1990 and 1995 [23. 

From the BCI data file it results that the studied area 
embodies 316 tree species, containing in total more than 
320,000 individuals The density of individuals for 

the 1995 census is illustrated on Fig. |4] with a gray-scale 



code (darker region corresponds to higher density). In 
the western part of the studied area lies a relatively large 
swampland (lighter region in the left side of Fig. [5j i.e. 
less trees). This area is ecologically quite different from 
the other parts (see for example the spatial distribution 
of the most abundance species: Hybanihus prunifolius 
(Fig. [5J. Therefore in determining the relevant macroeco- 
logical measures only trees located in the eastern square- 
like part of this region (500 x 500 meters) were considered. 




FIG. 5: Spatial distribution of the most abundant species: 
Hybanthus prunifolius. The square on the left shows the 
region in which the relevant macroecological measures were 
computed. 

The relative species abundance distribution for the 
census from 1990 is illustrated in Fig. [TJ The rank- 
abundance plot and the probability density curves sug- 
gest a scaling with an obvious cutoff. The scaling expo- 
nent is close to —1 (—0.98). Data for the other census 
years gives very similar results. 

For the BCI dataset the species-area scaling is not ob- 
vious at all, since the considered territory is quite small. 
As it is visible on Fig.[2]the power-law trend appears only 
in the limit of larger areas. The emergence of scaling be- 
comes more visible if one plots the local exponent (local 
slope) of the curves plotted in Fig. [2 Considering the 
census year 1995, on Fig. [6] we plot this variation. The 
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plot suggests in the limit of larger territories (area larger 
than lha) a stable power-law trend characterized by an 
exponent around 0.16 . As Fig. [2] already suggests, simi- 
lar results can be obtained for all census years, the stable 
scaling exponents varying between: 0.148 and 0.162. The 
obtained exponents are in agreement with species-area 
scaling results known in other macroecological systems 
0- 
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FIG. 6: Local slope (scaling exponent) of the curve plotted in 
Fig.[2]for year 1995. 

The C(r) auto-correlation function, characterizing the 
spatial distribution of the 10 most abundant species was 
also calculated. For all of them, a power-law type de- 
crease has been observed. As an example, on Fig. [3] 
the auto-correlation function of the species that is most 
abundant among trees with a stem thicker than 10 cm, 
Trichilia tuberculata, is visible (the spatial distribution 
of this species was illustrated on Fig. [5]). In Fig. [3K we 
plotted the results obtained on all trees belonging to this 
species and on Fig. [3b we plotted the results obtained 
for trees with diameter larger than 10 cm. Both picture 
suggest a power-law like decrease for C(r). 



V. A SIMPLE SPATIALLY EXPLICIT MODEL 

There are numerous theoretical modeling efforts for un- 
derstanding the relevant macroecological measures like 
the RSA distribution S EE EI EE El EI EE EE EE 
El, El, EE EE EE El El H- Most of these models are 
however mean-field like approaches in the sense that the 
spatiality of the individuals is lost, and the only relevant 
quantity that characterizes the system is the number of 
individuals in each species. These models are usually 
successful for reproducing the right RSA curves. Our 
goa l here is to follow the idea proposed by Zillio et al. 
[25l | to go beyond the mean-field approach and to study 
a simple but still realistic, spatially explicit model. The 
tropical tree system, where the coordinates of individuals 
are fixed, is especially suited for such an approach. 



The model considered by us is discrete both in time 
and space and it is inspired by the classical voter model 
[28| . Individuals are placed on a predefined lattice. Here, 
for the sake of simplicity a simple square lattice was con- 
sidered. Each individual belongs to a given species and 
the state variable characterizing each lattice site codes 
these species. It is assumed that the lattice with sizes 
L x L is always completely filled up, corresponding to 
a constant ecological saturation of the territory. Thus 
the total number of individuals in the system is always 
N = L x L. The number of possible states (species) is 
W, (W 3> 1). The dynamics of the system is governed 
by two adjustable parameters: p and g(0<p,g<l). 

In the beginning a randomly chosen Potts state 
(species) is assigned to each lattice site (individual). 
Starting from this random initial condition at each time 
moment a site is randomly chosen and its state is recon- 
sidered. This process models the disappearance of one 
tree and birth of a new individual in the opened niche. 
The p and q parameters will govern the selection of the 
new species that occupies this free position. 

• With probability 1 — q — p the state of the chosen 
site is changed to one of its 8 closest neighbors. 
This is the most probable process, since both p 
and q are much smaller than one. It models the 
quite common phenomenon, in which a seed from 
a neighboring plant will successfully develop into a 
tree. 

• With probability q a randomly selected species is 
assigned to the chosen site. This species is selected 
with a uniform probability from the ensemble of 
possible species. If W is big enough and q much 
smaller than one, this process is suitable for mod- 
eling speciation, or immigration in the considered 
territory. 

• With probability p a species already existing in the 
lattice is assigned to the chosen site. This species is 
selected by randomly choosing a site from the other 
N — 1 available ones and identifying the species at 
that site. The process models diffusion of seeds on 
the considered territory, allowing seeds originating 
from far away individuals to reach the considered 
location. 

In order to be in agreement with reality, the param- 
eters p and q are selected as: < q -C p <C 1. For 
p = and W ^> N we regain the analytically tractable 
model considered in [25j |. Due to the combined short and 
long-range interactions the analytical treatment of the 
model seems quite difficult. The model is studied thus 
by Monte Carlo (MC) simulations. In order to minimize 
the effect of boundaries, periodic boundary conditions 
are used. One MC step is defined as N updates of lattice 
sites. Simulations performed on several systems proved 
that in order to get realistic results for the final, statisti- 
cally stable state both q ^ and p ^ is necessary. For 
<7 = finally one of the species will colonize the whole 
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territory. On the other hand, for p = one would get 
unrealistically compact islands of species, in obvious con- 
tradiction with the fractal-like intermixing observed from 
the BCI dataset (see e.g. Fig. [7J. 
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FIG. 7: Spatial distribution of the abundant species. Fig- 
ure a. shows simulation results with p — 0, and figure b. is for 
p=0.03. Figure c. shows the spatial distribution of the abun- 
dant Trichilia tuberculata species in the BCI-1995 census. The 
other simulation parameters are L = 500, q = 1.5 x 10 . 



VI. MC SIMULATION RESULTS 

For q ^ and q <C 1 the model converges to a statisti- 
cally stable state with several coexisting species. As an 
example, considering the p =0 simple case and choosing 
L = 500, q = 0.0001 on Fig. we illustrate the time 
evolution of the total species number (W s ) in the sys- 
tem (thick continuous line). After 5000 MC steps the 
W s species number reaches a stable limit and fluctuates 
around W s « 400. On Fig. [8] we also show the time evo- 
lution of the population (JVj) for two selected species, one 
which appears at a later time moment and prevails, and 
one that appears quicker but gets extinct during the sim- 
ulation (thin dashed and continuous lines, respectively). 
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FIG. 8: Time evolution of the total number of species W s 
(thick continuous line) and the population of two selected 
species (continuous and dashed thin line). Simulations on a 
500 x 500 square lattice, q = 0.0001, p = 0. 

In order to get a realistic spatial distribution for the 
abundant species, p / is also necessary. For p = 
all species will be distributed in island-like structures 
(Fig.[7h.). Experimental results (Fig. Eh) show however a 



more homogeneous spatial distribution for all abundant 
species. By choosing a0<p€l value, one will get 
results in better agreement with the experimental ones 
(Fig. Eh). 

Several values for the p and q parameters (0 < p, q <C 
1) are appropriate for generating RSA and species-area 
curves in reasonable agreement with the experimental re- 
sults. Our goal here is to find those parameter values 
that will reproduce all the measured and previously dis- 
cussed features of the BCI dataset. First, we are looking 
for the p and q parameters that will lead to a statisti- 
cally stable species number — around 300 — in the simu- 
lated area. Second, it is also desirable to get the relative 
size of the most abundant species (number of individuals 
in the most abundant species divided by N) of the same 
magnitude as in the experimental data. For the BCI 
dataset (census years 1982, 1985, 1990 and 1995) this 
relative size is: 0.169, 0.17, 0.165 and 0.157, respectively. 
Third, we are looking for a species-area scaling exponent 
in the 0.14-0.17 range (Fig. HJ and an RSA probability 
density that shows a scaling with a —1 exponent and an 
obvious cutoff for abundant species (Fig. [lb) . 

An extensive search by MC simulations in the p-q pa- 
rameter space suggested the p — 0.03 and q — 1.5 x 10 -5 
values. For these parameter values we got a stable species 
number around 300 and the relative size of the most 
abundant species 0.14. The species-area curve (Fig. [9]) 
suggests a scaling exponent of 0.156. 
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FIG. 9: Simulation results for the species-area curve (log-log 
scale) after 10000 MC steps. The power-law fit in the limit of 
large areas is drawn with a dashed line. The scaling exponent 
is 0.156. The unit for area is chosen as 100 x 100 lattice sites. 
Simulation parameters: L — 500, p = 0.03 and q = 1.5 x 10~ J . 

Moreover, finite-size effects in the species-area curve 
shape are also very similar with the ones experimentally 
obtained. Plotting the local slope of the curve in Fig. [9[ 
the trend presented in Fig.QJJis obtained. This trend can 
be compared with the one suggested by the BCI dataset 
(Fig. [6] for the 1995 census data). 

The shape of the RSA curve is also in good agreement 
with the BCI results. Both the Preston plot (Fig. Hlk ) 
and the probability density function (Fig. [TTb) resembles 
the experimentally observed ones (Fig. [J). The p(s) prob- 
ability density has a scaling in the limit of rare species 
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FIG. 10: Local slope (scaling exponent) of the curve plotted 
in Fig. [9] 



characterized by an exponent —0.97, close to the experi- 
mentally observed —0.98 value (Fig.QJ;). 

Unfortunately, this spatially explicit model does not re- 
produce the observed power-law type decay (Fig.[3| of the 
auto-correlation function. The model suggest an expo- 
nential decay for the abundant species for all reasonable 
p and q parameter values. In particular, for the above 
discussed optimal p and q values, again an exponential 
decay is observed for the abundant species (Fig. [ 
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VII. DISCUSSION AND CONCLUSIONS 
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A simple spatially explicit model, which was first con- 
sidered in [23], was extended here for explaining the sta- 
tistical properties of individuals and species in a tropical 
tree community. The model is defined on a square lattice 
and it is inspired by the classical voter model [28]. It de- 
scribes in a realistic manner the spatiality of individuals, 
birth and death processes, diffusion of seeds and specia- 
tion events. Within the model the dynamics of species 
is governed by two adjustable parameters, p and q. The 
parameter q defines the probability of speciation, or im- 
migration events. The parameter p describes the proba- 
bility for the global seed diffusions within the considered 
territory. 

The system was studied by large scale Monte Carlo 
simulations and the p and q parameters were adjusted for 
reproducing optimally all experimentally studied macroe- 
cological measures. For reasonable values of p and q the 
model is successful in explaining the experimentally de- 
tected species number, relative size of the most abun- 
dant species, the good species-area scaling and the shape 
of the RSA curves. The model also generates a visually 
good spatial distribution of individuals within a species. 
The best model parameters were found as p = 0.03 and 
q = 1.5 x 10~ 5 . Beside all these successes, the model fails 
to reproduce the experimentally observed power-law de- 
cay for the auto-correlation function of the individuals 



FIG. 11: Simulation results for the RSA curves after 10000 
MC steps. Figure a. shows the classical Preston plot, figure b. 
the rank-abundance curve and figure c. the p(s) probability- 
density function. The dashed line suggests a power-law with 
exponent — 1. Simulation parameters: L = 500, p = 0.03 and 
g=1.5xl0" 5 . 



within a species as a function of separation distance. As 
one would naturally expect from a lattice models where 
the local interactions are dominating (p,q <C 1), simula- 
tions lead to an exponential decay. This result suggests 
that in reality long-range interactions might be more im- 
portant than it was considered in this simple approach. 

As final conclusion, it can be stated that in spite of its 
simplicity, the present model is successful in explaining 
the majority of the experimentally available macroeco- 
logical measures. More complex models are needed how- 
ever in order to fully understanding all aspects of spatial 
distribution of species in the tropical tree communities. 
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